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We discuss the generation and statistics of the density fluctuations in highly compressible poly- 
tropic turbulence, based on a simple model and one-dimensional numerical simulations. Observing 
that density structures tend to form in a hierarchical manner, we assume that density fluctuations 
follow a random multiplicative process. When the polytropic exponent 7 is equal to unity, the local 
Mach number is independent of the density, and our assumption leads us to expect that the prob- 
ability density function (PDF) of the density field is a lognormal. This isothermal case is found to 
be singular, with a dispersion a 2 which scales like the square turbulent Mach number M 2 , where 
s = In p and p is the fluid density. This leads to much higher fluctuations than those due to shock 
jump relations. 

Extrapolating the model to the case 7 / 1, we find that, as the Mach number becomes large, 
the density PDF is expected to asymptotically approach a power-law regime, at high densities when 
7 < 1, and at low densities when 7 > 1. This effect can be traced back to the fact that the 
pressure term in the momentum equation varies exponentially with s, thus opposing the growth of 
fluctuations on one side of the PDF, while being negligible on the other side. This also causes the 
dispersion a 2 to grow more slowly than M 2 when 7 / 1. In view of these results, we suggest that 
Burgers flow is a singular case not approached by the high-Af limit, with a PDF that develops power 
laws on both sides. 
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I. INTRODUCTION 

The formation of density structures by the velocity field of highly compressible turbulence is of great interest in 
astrophysics. The determination of their typical amplitude, size and volume filling factor poses significant difficulties 
since it requires a knowledge of the full statistics. As a first step, we shall concentrate in this paper on one-point 
statistics and more specifically on the probability density function (PDF) of the density fluctuations in one-dimensional 
(ID) turbulent flows. 

It is well known that the density jump in a shock depends directly on the cooling ability of the fluid. Thus, for an 
adiabatic flow the maximum density jump is 4, for an isothermal flow it is ~ M% p|, and for nearly isobaric flows it 
is ~ e M » |Q, where M a is the Mach number ahead of the shock. The net cooling ability of a flow can be conveniently 
parameterized by the polytropic exponent 7, so that the thermal pressure P is given by P = Kp 1 , where p is the fluid 
density ||. Isothermal flows have 7 = 1, and isobaric flows have 7 = 0. Note that 7 < corresponds to the isobaric 
mode of the thermal instability (see, e.g., Q). Thus, in general, the amplitude of the turbulent density fluctuations 
will be a function of 7. 

Previous work with isothermal flows had suggested that the PDF is log-normal (111, while for Burgers flows a 
power-law PDF has been reported 0. More recently, evidence that flows with effective polytropic indices < 7 < 1 
also develop power-law tails at high densities has been presented || . In order to resolve this discrepancy, we present a 
series of ID numerical simulations of polytropic gas turbulence with random forcing, in which the polytropic exponent 
7 parameterizes the compressibility of the flow. We have chosen to use ID simulations in order to perform a large 
number of experiments at a sufficiently high resolution, integrated over very long time intervals, allowing us to collect 
large statistical samples. 

The simulations have three governing parameters: the polytropic index 7, the Mach number M, and the Reynolds 
number R. We keep the Reynolds number fixed, and explore the effects of varying 7 and M on the resulting density 
PDF. We find that varying these two parameters is not equivalent. Variation of 7 induces a clear qualitative variation 
of the density PDF, which, at large Mach number, displays a power-law tail at high densities for < 7 < 1, becomes 
log-normal at 7 = 1, and develops a power-law tail at low densities for 7 > 1. This suggests a symmetry about 
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the case 7=1, which we also explore. Variation of the Mach number, on the other hand, only appears to induce a 
quantitative change, in such a way that increasing M augments the width of the PDF. 

T he p lan of the paper is as follows. In sec. O we describe the equations solved and the numerical method. In 
sec. |n| we describe the statistics of the various fields, in terms of their PDFs, together with a tentative model and 
a discussion of the Burgers case. Section |^ is devoted to a discussion on the choice of the forcing, together with a 
summary of our results. 



II. EQUATIONS AND NUMERICAL METHOD 

We choose to concentrate on one-dimensional forced polytropic gas dynamics, governed by the following non- 
dimensionalizcd equations 

d t u + ud x u = ]—r^— + ^-d xx u + a (1) 

"fM A p R 

d t p + d x (pu)=0 (2) 

where u is the velocity of the fluid in units of U, p the density in units of po, 7 the polytropic index and M the Mach 
number of the unit velocity U at the unit density p Q . The equations are driven by an acceleration a with zero mean. 
The Reynolds number is R = where L is the size of the domain and v the kinematic viscosity chosen constant to 
ensure the conservation of the mean velocity (u) = -r J udx. The viscous term is kept as small as possible and is only 
here to prevent numerical blow-up. Note that the "correct" form of the viscous term is obtained after replacing v by 
the ratio pj ' p, where the dynamical viscosity p is usually considered independent of the density. The equations then 
conserve the momentum J pudx if the acceleration a in eq. (|l|) is also replaced by the ratio of a force / to the density 
p. The dynamics that results in this case is very different due to the dependence of the driving term with respect to 
the density, as discussed in the last Section. 

For large Mach number simulations, it was found necessary to smooth density gradients, using a mass diffusion 
term of the form p r d xx p in the right-hand side of eq. (^|). Total mass is still conserved in the presence of this term, 
and if p r is taken sufficiently small, it has been tested that it does not affect the dynamics in a way that could modify 
our conclusions. 

We also found it convenient to solve eqs. (l])-(Q) using the variable s = In p. The numerical code uses a standard 
pseudo-spectral method with periodic boundary conditions. Time advance is performed using a Crank-Nicholson 
scheme for the linear terms and an Adams-Bashforth scheme for the nonlinear ones. For all the runs presented in this 
paper, the kinematic viscosity has been fixed to v = 3 x 1CP 3 . For runs with M > 3, we have /i r = 5x 10~ 4 . 

The acceleration a is prescribed in Fourier space. Its spectrum has a constant amplitude (equal to 0.6) on wavenum- 
bers 1 < k < 19 and phases chosen randomly with a correlation time t COI = 0.003. Resolution ranges from N — 3072 
to N — 6144 grid points for the runs with M > 6. 

We perform one point statistics of the simulations, both for the density and the velocity derivative, keeping the 
forcing and the viscosity constant. All simulations start with zero initial velocity and constant density. 

In order to obtain reasonably sampled histograms of the one-dimensional fields, which contain only N spatial data 
points, we sum the histograms over time, sampling at intervals of 0.1 time units, integrating over a total of 150 time 
units. However, we have found that, since the simulations start with uniform density, the first several samples must 
be discarded, since they bias the density histogram near p = 1. We typically skip the first 20 temporal samples (2 
time units). The PDFs thus computed contain roughly 4 million data points. Note that longer integration times 
are needed at larger Mach number in order to reach a statistically relevant sample, the sound crossing time of the 
integration domain being larger as M increases. 



III. A MODEL FOR THE DENSITY PDF 

A. Properties of the governing equations 

Before describing our model for the density PDF, it is instructive to rewrite the governing equations in the inviscid, 
unforced case, using the variable v = (1 — 7) hip when 7^1 and s — hip when 7 = 1. We get, for 7 =/= 1 

— = 1 9 c-v (3) 
Dt (l-^M^dx { > 

Dv d 
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and for 7=1 

Du 1 d 

Dt~ APdx S {) 

Di = -d~ X U (6) 

where -St stands for the convective derivative. The variable v is, up to an additive constant, the logarithm of the 
square of the sound speed, and when 7 = 1, becomes identically zero. 

These equations can be rewritten in Ricmann invariant form. For 7^1, they read 

[d t + (u±c)d x ](u±-r^ v -) = 0, (7) 

7 — 1 

where c = /M is the sound speed, while in the singular case 7 = 1 these equations become 

[ dt + (u±-L)d x ](u± 1 ^)=0. (8) 
A number of interesting remarks can be made on the previous equations. 

(i) When 7=1, eqs. (^)-@ are invariant upon the change s — > s + b, where b is an arbitrary constant. Indeed, the 
sound speed does not depend on the local density of the fluid. 

(ii) In the general case, if we substitute 7 by 2 — 7 and p by we observe that the Riemann invariants z^ = 
u ± (~!x\ & re exchanged, while their speeds u ± c remain unchanged. We shall now explore the implications of this 
remark on the statistics of the density fluctuations in the weakly compressible regime. For small values of the Mach 
number^ a reductive perturbation expansion can be performed on the viscous equations and it has been shown Q (see 
also [[10|) that one-dimensional compressible turbulence reduces essentially to the superposition of the solutions of 
two Burgers equations describing nonlinear wave propagation in opposite directions. More precisely (considering eqs. 
( |l|j^ ) with M = 1), if we denote by p' and vl the perturbations of the basic state (p = 1, u = 0), Tokunaga obtained 

^ = -^(fi(ei,r)-F 2 (fc,r)) (9) 
7+1 

u' = -^-r(F 1 (Z U T) + F 2 (b,T)) (10) 
7+1 

where e is the order of magnitude of the nonlinear waves. The new coordinates and r arc defined by 

& = e[x - nt- <pi(x,t)] (11) 
r = e 2 t (12) 



where n = 1 and T2 — — 1 and the phase functions obey 

13-7 
21+^ 
13-7 



21 + 7 



fi(£,r)de + 03, (14) 



with 9i arbitrary constants determined by the initial conditions. Finally the functions Fi (simply related to the 
Riemann invariants z^) satisfy the Burgers equations 

d T F i +F i dt t F i = ~dte i F i . (15) 

The fields Fi evolve almost independently, with the same dynamical equation, except for phase shifts, a higher order 
effect most important during collisions of shock waves. Given some initial conditions for p' and u', the substitution 
p — > 1/p (or p' — > —ff), and 7 — > 2 — 7 leads to the replacement of F\ and F 2 by Fa(3 — 7)/(l+7) and jPi(3 — 7)/(1+t) 
respectively. For a vanishingly small viscosity v, the rescaling of the amplitudes F\ and F 2 can be absorbed in a 
rescaling of the variables Except for this stretching of the space and time variables, this substitution will thus lead 
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to the same fluctuations occurring at different locations. As a consequence, we can expect that the probability density 
functions of the cases 7 and 2 — 7 for small values of M will be closely related after the change p — > 1/p. The case 
of higher Mach numbers is more delicate due to the additional problem of mass conservation, rendering impossible a 
symmetry between p and 1/p. This question is addressed below. 

(iii) The substitution 7^2 — 7 can also be examined at the level of eqs. (||)-([|). Its effect is simply to change the 
sign of the right-hand sides. For 7 < 1, eq. (|j) shows that positive values of v (in this case associated with density 
peaks) are mostly created by shocks (associated with negative velocity gradients). Looking at eq. (||), we see that 
as v increases, the pressure term becomes exponentially small and thus cannot prevent the formation of very strong 
peaks. Negative values of v (here associated with density voids) are created by expansion waves, but in that case 
the pressure increases exponentially with decreasing values of v leading to a rapid saturation of this process. As a 
consequence, we expect that for 7 < 1 the PDF of v will be significantly more populated at positive rather than at 
negative values. For 7 larger than unity the PDF of v will be similar, the formation of positive values of v (now 
associated with density voids) being still unhindered by the pressure. It results that the PDF of s = In p for 7 > 1 
will appear similar to that for 7 < 1 after we change s — > — s. 

(iv) When 7=1, the behavior is very different since the acceleration due to the pressure term is simply proportional 
to —d x s and thus never becomes negligible. We expect a symmetry in the PDF of s, positive and negative values of 
s being equally created, by shocks and expansion waves respectively. 

(v) Finally, it is also useful to discuss the shock jump relations for the polytropic equations. Denoting by X the 
ratio of the post-shock to the pre-shock density, we have (see j|] ) 

X 1+1 - (1 + 1 m 2 )X + 1 m 2 = (16) 

where m is the upstream Mach number in the reference frame of the shock. This equation shows that for 7=1, 
X = m and that the jump X increases more slowly than m 2 for 7 > 1, while it increases faster than m 2 for 7 < 1 
with, as 7 — > 0, X ~ e m . For weak shocks, we get X r* 1 + (m 2 — In this case, the shock velocity being 

close to the sound speed, we can write m — 1 + -j-, where u is the velocity in the simulation frame, leading to 
X = l + ^- = l + -j^— ■ We thus get ^ ~ m s j^, with m s = u/c denoting the Mach number in the simulation 
frame. 

Note that typical pressure fluctuations created by almost incompressible turbulence scale like M 2 where the tur- 
bulent Mach number is defined as M — u rms /c (here u rms is mostly made of solenoidal motions unlike in our ID 
simulations where it stands for purely compressible modes). This scaling corresponds to a balance between the pres- 
sure gradient and the nonlinear term. If entropy fluctuations are not allowed (like with a polytropic state law), the 
resulting density fluctuations also have to scale as M 2 . In thermally forced turbulence however, a Boussinesq-like 
balance obtains between temperature and density fluctuations, maintaining pressure fluctuations of order M 2 , while 
allowing for much larger values of density and temperature fluctuations JTl| . 

In weakly nonlinear acoustics, the pressure term is balanced by the velocity time derivative and we recover the 
scaling — ~ M obtained for weak individual shocks. 

B. The case 7=1 

The main idea of our model is that density fluctuations are built up in a hierarchical process [||. After a shock 
(respectively an expansion wave) passes through a given region of mean density po, the density reaches a new value 
Pi, larger (respectively smaller) than p a . In this region new fluctuations can be created, changing the local value p\ 
to P2 and so on. Of course the dynamical equations constrain this process. For example, due to mass conservation, 
arbitrarily high values of the density can only be reached in very localized and thin peaks. We thus expect this 
hierarchical process to saturate at some value s+ > 0. A similar saturation should occur for low densities at some 
value s_ < 0, with probably |s_| > |s+|, due to the fact that larger voids can be created without violating the mass 
conservation constraint. 

The build up of these density fluctuations is a random multiplicative process which, at the level of the variable s, 
becomes additive. The random variable s is thus the sum of individual random variables (the density fluctuations), 
each having the same probability distribution. The latter fact follows from the invariance of the equations at 7 = 1 
under the change s — > s + s , which furthermore implies that each individual jump has the same average magnitude, 
related to the Mach number of the flow but independent of the local density. The sum of identical random processes is 
known to have a Gaussian distribution, due to the Central Limit Theorem, whatever the distribution of the individual 
processes. The PDF of s is thus expected to follow a normal distribution. 

The variance of the random variable s can be estimated using the size of the typical fluctuations associated both to 
shocks and expansion waves. The case of shock waves has been discussed above. At small values of M, Pi+i/ Pi ~ 1+M 
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(with the Mach number M = u rms /c ) so that Ss = ln(p i+1 /pi) = ln(l + 5p/p) w ln(l + M) - M. At high Mach 
numbers, the individual jumps obey As ~ InM 2 . For expansion waves, the balance of the time derivative of s and of 
the positive velocity gradient in eq. (||), gives s ~ M, regardless of the value of M since the term us x is smaller than 
the term s t if the density decreases uniformly in space (note that the decrease of p is exponential in time). Thus, 
with this mechanism, the density decreases in the center of the expansion waves while it increases on the edges, until 
pressure blocks the process. In the case 7 = 1 pressure acts symmetrically in s and we thus get positive and negative 
fluctuations which are of the same order of magnitude, and much larger than those due to shocks. We thus expect 
that a s ~ M for a large range of values of the Mach number. 

From the previous dicussion we can expect the PDF of the variable s to be given by 

P(s)ds = _L=exp(-(^!V (17) 

with <r 2 = /3M 2 , and (3 a proportionality constant. The maximum of this distribution s is simply related to a s due 
to the constraint of mass conservation. Writing (p) = J_°° e s P(s)ds — 1 we find s — — icr 2 (see below). Note that 
the PDF of p is related to that of s by P p (p) = P(ln p) /p. 

The predictions of this model can be tested against results from numerical simulations. Figure [l] (top panel) shows 
a plot of log(<7 s ) vs. log(M) obtained by combining data from several simulations with M — 0.5, 1, 2, 3, 4.5, 6 and 10. 
These data were obtained by computing M and a s for the accumulated density and velocity fields over 100 subsequent 
outputs of the simulations (spanning a total duration of 10 time units) for each point in Fig. ^. This plot shows that 
<7 2 ~ f3M 2 , with P w 1, with a very good accuracy, up to the highest Mach numbers reached in our simulations. 
On the other hand, we see in the bottom panel of the same figure, which displays log<7 p vs. logM, that the density 
standard deviation also scales like M for small values of M , while for M > 0.5 the points curve up, a reflection of the 
relation er 2 = — 1 between the two variances when p obeys a log-normal distribution. The relation s a — — i<r 2 is 
also well verified numerically as can be seen form Fig. g. 

We now display in Fig. || the logarithm of the s-histograms for three runs with 7 = 1 and M — 0.5, 2, and 6. Fits 
with parabolas are shown in dashed lines and show that, to a very good approximation, the PDFs of the density are 
in all three cases log-normals. An estimation of the widths and maxima of these distributions also shows a very good 
agreement with the predictions a s ~ M and s a = — 0.5a 2 . 

The distribution of the velocity derivative u x is shown in Fig. [|for 7 = 1 and M = 6. This distribution is found to 
be almost independent of the Mach number. It presents a long exponential tail for negative values of u x and a strong 
drop off for large values, analogous to the one found in the Burgers case 0. 



C. The case 7 / 1 



The difference between the case 7 = 1 and the cases 7 ^ 1 lies in the behavior of the pressure term as a function 
of the local density of the fluid, an effect which is most visible after comparing eq. (^) with eq. (||). With the 

density-dependent rescaling M — > M(s;~f) — MeT 2 *, the two equations identify, which only means that we expect 
the small-fluctuation behavior of the case 7 ^ 1 to be identical to that of the case 7 — 1 in regions where the local 
logarithm of the density is close to s, when M(s; 7) is substituted for M. 

The argument at the origin of the PDF of s in the isothermal case is based on the fact that the local Mach number 
of the flow is independent of the local density. When 7^1 this property is violated and there is no reason to expect 
a log-normal PDF for the density. We nevertheless propose a heuristic model, reproducing most of the features of the 
PDF's obtained in our simulations, which consists in taking the same functional form of the PDF as in the isothermal 
case, but replacing M by M(s;~f), where M(s;j) now stands for the "effective" r.m.s. Mach number at the value 
s. This "effective" r.m.s. Mach number is defined as M(s\^f) = u rms /c(s), the local turbulent Mach number, when 
s_ < s < s + , and by the constant M(s + -j) (respectively M(s_;7) ) for s > s + (respectively s < s_). These cutoffs, 
which, as we shall see, are necessary for convergence, are also physically meaningful, since the probability of new 
fluctuations arising within previous peaks or voids decreases as the amplitudes of the latter become larger because the 
fraction of space they occupy decreases. The fact that the cutoff occurs at larger values of |s| for s < than for s > 
is due to the larger filling factors of low density regions (see Fig. ^a and Fig. |^b for comparison). A numerical check 
of this saturation property is possible if one computes the scatter plot of the standard deviation for s vs. the mean 
value of s in subregions of the integration domain for each snapshot. Figure ^ shows these plots for M = 6, 7 = 0.5 
and 7 = 1.5 in subregions of length N/3. A clear trend is visible, indicating the change of the local Mach number 
with the local mean density. Moreover, we clearly see that the saturation level for s < at 7 = 1.5 occurs at a much 
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higher value of the Mach number than for s > 0, 7 = 0.5. Plots of a s and a p vs. M for 7 = 0.5 and 7 = 1.5 are also 
presented in Figs. |t] and |[ They show that a s increases more slowly than linearly with M for high Mach numbers. 
This results from the asymmetry in the fluctuations of s for 7^1. While for 7 = 1 the typical excursions of s are 
of the order of M both for positive and negative values of s, when 7 > 1 for example, pressure blocks the negative 
fluctuations of s while still allowing for fluctuations of order M on the positive side. The resulting variance a s is thus 
expected to be smaller than M. The same argument applies for 7 > 1 but then fluctuations are of smaller magnitude 
when s > 0. Looking at the plot of a p wc see opposite trends for 7 > 1 and 7 < 1. First, note that we do not expect 
the specific relation mentioned above between the variances of s and p, since the distribution of s is not Gaussian. 
Second, this trend is easily interpreted if we recall that for 7 < f the density fluctuations are in high peaks, while for 
7 > 1 they consists of large voids. In the former case the variance of p can increase greatly when M is large, while in 
the latter case, the voids do not contribute much in the variance, leading to a slower increase of a p with M. 
The PDF will thus read 



s 
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2M 2 ( S ; 7 ) 



P(s;~/)ds — C(-/) exp — - — a(j)s ds = (7(7) exp 2 A/ 2 ct{"f)s ds (18) 



2 p (7-l)s 



-s e 



where C(j) is a normalizing constant such that J^°° P(s;^y)ds = 1. The parameter 0(7) is again determined by the 

constraint of mass conservation stating that the mean value of the density should be 1 : f_°° e s P(s;'f)ds = 1. Note 
that in the absence of cutoffs, the convergence of the integrals would require a > 1 for 7 < 1 and a < for 7 > 1. This 
functional form of the PDF immediately allows to make a few predictions. For 7 < 1, Af(s;7) grows exponentially 
with s for s- < s < s+ and as a consequence the PDF for < s < Af (s; 7) is dominated by the power-law (in p) 
behavior P(s;j) ~ e^ a ^ s , while the Gaussian-like decay will again dominate for s > M(s;j). For s < 0, the local 
turbulent Mach number decreases as s decreases and we expect a drop off of the PDF more rapid than when 7=1. 
The behavior is exactly opposite when 7 > 1. This prediction can be verified by looking at Fig. displaying the PDF 
of s for 7 = 0.3 and 7 = 1.7 at M — 3. 

It is now interesting to relate the PDF for a certain value of 7 to that obtained for 2 — 7. Writing the condition 
(p) = 1, we get 



exp 



s 2 \ f 00 ( s 2 \ 

— - h (1 - aM)s )ds = / exp aMs)ds (19) 

2Af 2 ( S ; 7 ) UJJ ) J^oo V 2M 2 (s; 7 ) K1> ) V ' 

while the same condition for 2 — 7 reads, after making the substitution s — ► — s in the integrals 

s 2 r°° s 2 

exp( R (l-a(2-^))s)ds= exp( s h a(2 ~ j)s)ds. (20) 

, V 2M 2 (- S ;2- 7 ) V , " ' Loo V 2Af 2 (-s;2- 7 ) "' 

For s_ < s < s + , the functions Af(s;7) and Af(— s; 2 — 7) are identical. If the cutoffs s+ and s_ occur at large 
enough values, i.e. when the local Mach number is either very large or very small, the contributions in the integrals 



of the two terms involving these two quantities will be very close and, by inspection of eqs. (19) and (EOl) we get 



a(2- 7 ) = l-a( 7 ). (21) 

This relation is exact when 7 = 1 since M(s; 1) = M is independent of s, allowing to recover the result a(l) = i. 
Note also that for large enough M, a case where eq. (^lj) holds, the symmetry s — > — s is not possible but must 
include a translation in the s domain to account for mass conservation. 



Relation (21) is verified numerically with a reasonable precision for the runs at the highest Mach numbers. For 
example, when M = 6, the slope of the power law is —1.2 (i.e. a — 1.2) for 7 = 0.5 while we have a = —0.28 for 
7 = 1.5 (see Fig. |l^). For smaller values of M, the absolute values of the slopes are closer to each other, a feature 
due to the different cutoffs for negative and positive values of s (see Fig. || for M = 3 and 7 = 0.3 and 0.7). Note 
that the shape of the PDF for M = 6, 7 = 1.5 presents a steeper slope for values of s slightly smaller than that of 
the maximum. This feature can also be reproduced with this simple model, as can be seen on Fig. [ll] which displays 
the PDF obtained from eq. (jT|) for a = 0.28, 7 = 1.5, M (0) = 1.2 and values of the cutoffs at M = 10 for s < and 
M = 0.1 for s > 0. 



D. The case 7 = 0, i.e. Burgers' equation 

An interesting problem concerns the high Mach number limit. It is often suggested that when M is very large, the 
dynamics of compressible flows should be analogous to that prescribed by the Burgers equation. While this may be 
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true for the velocity field, our results prove that it cannot be the case for the density. Indeed, we find that whatever 
the value of 7 ^ and of the Mach number, there is always a range of densities for which the pressure cannot be 
neglected. For that range of densities the PDF has no power- law tail but presents a more rapid drop-off. For 7=1, 
it turns out that the pressure is never negligible. Extrapolating our results, we thus predict that for the Burgers' 
case there should be power law tails both for low and high densities. We thus performed a simulation of the Burgers 
equation (coupled with eq. (|j) for the density) with the same parameters as for the previous runs and with N = 6144. 
The resulting PDF is presented in Fig. This plot shows that indeed the PDF is almost flat for s < 0, while there 
is also a power law for s > 0, with a negative slope of roughly 0.5. The cutoff for large densities is due to the viscous 
terms, which give a minimum scale for the width of the shocks, and thus a maximum value for the density peaks. In 
the physical domain, we observe the creation of voids (s reaching a value of —85 at t = 64) which occupy most of 
the domain, together with very high peaks (s rs 6.5). The number of peaks decreases during the simulation while the 
density in the voids decreases exponentially in time. The forcing is unable to break the peaks because it acts at large 
scales, while the the density fluctuations become as narrow as allowed by viscosity. This PDF has to be contrasted 
with the one obtained in Q, for which the Reynolds number was low and the simulation decaying. In that case the 
power law at high densities was obtained but the PDF presents a sharp drop off for low densities. Two-dimensional 
decaying simulations of the Burgers equation are also presented in Q for moderate Reynolds numbers. The plateau 
of the PDF at low values of the density is also obtained, while the power law for s > is not as clear. Burgers 
simulations for the decaying infinite Reynolds number case are presented in |l2| . In that case the PDF is calculated 
for the cumulated mass function and not for the density, which is not defined after the first shock formation. A power 
law is found, which extends to s — —00 and connects to an exponential decay for s —> +00. Note that an exponential 
PDF for the density was predicted in jl3| on the basis of a model which treats shocks as completely inelastic particles. 
We can thus conclude this section by saying that the Burgers case is truly a singular limit, which cannot be reached 
as the high Mach number limit of a polytropic gas, with 7^0. 

IV. DISCUSSION 
A. Effects of the forcing 

The study presented in this paper has been performed for a single choice of the forcing and of the Reynolds number. 
While the variation with the latter parameter can be trivially extrapolated, we cannot a priori be sure that our results 
are independent of the type of forcing. We have performed decay runs and observed that the behavior of a s vs. M 
is still the same as in the forced case. The PDFs however cannot be computed on a single snapshot due to the poor 
statistics and cannot be integrated in time since the Mach number changes by roughly one or two orders of magnitude 
during the run. We have also performed a run at 7 — 1 with a forcing of the form f / p in eq. (Q). In that case 
the density PDF is not a lognormal anymore but presents a power law tail for low densities (not shown). This can 
be attributed to the fact that the flow is stirred more vigorously at low densities so that the effective Mach number 
indeed increases as p decreases. We nevertheless think that our results can be extrapolated to an unforced situation, 
at a given time, and possibly also to the multi-dimensional case. Note that the Mach numbers we have explored in 
this paper would correspond to even higher Mach numbers in the multi-dimensional case since in that case only a 
fraction of the total kinetic energy populates the compressible modes. 

B. Summary 

We have presented an investigation of the density PDFs of a randomly accelerated polytropic gas for different values 
of the polytropic index and of the Mach number. We have suggested a simple model in which the density field is 
everywhere constructed by a random succession of jumps j^]. When the flow is isothermal (7 = 1), the jumps are 
independent of the initial density, and have always the same probability distribution. Expressed with the variable 
s = In p the jumps are additive, and by the Central Limit Theorem are expected to have a Gaussian PDF, or a 
lognormal in p. 

An analysis of the expected s increments in the weak and strong shock cases, as well as those due to expansion 
waves, suggested that the variance a 2 should scale as the mean square turbulent Mach number M 2 . Moreover, because 
of mass conservation, the peak of the distribution s a is related to the variance by s a — — \<J 2 . These predictions were 
verified in ID simulations of compressible turbulence. Previous claims that it is the density variance a 2 that should 

scale as M 2 || might have been misled by lower effective Mach numbers than those achieved in the present simulations, 
in which all of the kinetic energy is in compressible modes thanks to the one-dimensionality. 
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When 7^1, the density jumps are not independent of the local density anymore, and the shape of the PDF should 
change. Observing that a renormalization of the Mach parameter (eq. (jl])) M — > M(s; 7) = Me~^ Ls restores the 
form of the equations for the case 7 = 1, we proposed the ansatz that the PDF may still be described by the same 
functional form as in the case 7=1, but substituting M by M(s;7). This prediction is confirmed by the numerical 
simulations, giving PDFs which are qualitatively in very good agreement with the model PDF, eq. The result 

is that the PDF asymptotically approaches a power law on the side where (7 — l)s < 0, while it decays faster than 
lognormally on the other side. 

Upon the replacements 7 — > (2 — 7) and p — > 1/p we find, using the condition of mass conservation, that the slope 
a of the power law for a given value of 7 is related to its value at 2 — 7 by eq. (|2l]) in the large Mach number limit. 
These results are also confirmed by the numerical simulations, which exhibit a power law at s > when 7 < 1 and at 
s < when 7 > 1, with slopes which are roughly related by eq. (|2l|), with better accuracy at large Mach numbers. 

Finally, on the basis of these results, we suggested that the Burgers case should develop a power law PDF at both 
large and small densities, since in this case there is no pressure on either side. This result was again confirmed by a 
simulation of a Burgers flow. 

We shall conclude this paper by pointing out that the non-uniqueness of the infinite Mach number limit might have 
important consequences for astrophysical applications, such as in cosmology. The so-called Zeldovich |l4j approxima- 
tion is indeed based on the Burgers' equation which, at the light of the present work, appears as a questionable model 
of highly compressible flows. This point will be addressed in future work. 
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Figure captions 

FIG. 1. (Top) Variance of s = hip vs. the mean square Mach number M 2 mE = M 2 for various simulations with 7 = 1 and 
M — 0.5, 1, 2, 3, 4.5, and 6. Every point in this plot gives the variance and M for sets of 100 subsequent outputs (10 time 
units) of any given simulation. The simulations were typically run for 150 time units. (Bottom) Variance of p vs. M. 



FIG. 2. Most probable value of s vs the variance of s, a 2 , for the runs in Fig. Ftl The data points are obtained as in Fig. [jl 
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FIG. 3. Probability density function (PDF) of s for three simulations with 7 = 1 and M = 0.5, 2 and 6. For clarity, these 
PDFs have been respectively displaced in the plot by —2, — 1, and units in the vertical axis. The shift of the peak towards 
more negative s values at larger M is real, due to the constraint of mass conservation. The dashed lines show the best fit with 
a lognormal to each PDF. 



FIG. 4. PDF of the velocity derivative for a run with 7 = 1 and M — 6. 



FIG. 5. a) (Top) Density field of a run with 7 = 0.5 and M = 10 at time t = 34.65 Note the very thin density peaks and 
the shallow density minima, b) (Bottom) Density field of a run with 7 = 1.5 and M — 6 at t = 50.5. Note that the density 
maxima are now much shorter, while the density minima (voids) become much deeper. They are also much wider than the 
peaks in the 7 = 0.5 case because of mass conservation. 

FIG. 6. Standard deviation of s vs. the mean value of s over subregions of size 1/3 of the integration domain for two runs 
with (top) 7 = 0.5 and (bottom) 7 = 1.5. Note the inverse trends between the two runs and the saturation of a s at large values 
of |(s)|, especially noticeable in the case 7 = 1.5. 



FIG. 7. Variance of s (top) and of p (bottom) vs. the mean square Mach number for 6 runs with 7 = 0.5 and M = 0.5, 2, 3, 
4.5, 6 and 10. Note that ct 2 increases more slowly than M r 2 m8 because only one side (s > 0) of the density PDF is unimpeded 
by the pressure. Instead, cr 2 increases more rapidly than M r 2 ms because such fluctuations in s imply very large fluctuations in 
P- 



FIG. 8. Variance of s (top) and of p (bottom) vs. the mean square Mach number for 6 runs with 7 = 1.5 and M = 0.5, 
2, 3, 4.5 and 6. Again (compare to Fig. [?]), <r 2 increases more slowly than M 2 ms because only one side s < of the PDF is 
unimpeded by the pressure. However, in this case also a 2 increases more slowly than M 2 ms , because the density fluctuations 
are bounded by zero, not being able to contribute much to the variance of p. 

FIG. 9. PDFs of s for two simulations with M — 3 and 7 = 0.3 (top) and 7 = 1.7 (bottom). For 7 = 0.3 the power-law 
regime appears at large densities, while for 7 = 1.7 it appears at small densities. 

FIG. 10. PDFs of s for two simulations with M = 6 and 7 = 0.5 (top) and 7 = 1.5 (bottom). Note that at this Mach 
number, the power-law regime for the 7 = 1.5 case appears removed from the peak of the distribution, mediated by a regime 
with a steeper slope. 



FIG. 11. The theoretical PDF given by eq. (j 
s < and M = 0.1 for s > 0. Compare to Fig. 1 



for a = 0.28, 7 = 1.5, M(0) = 1.2 and values of the cutoffs at M = 10 for 



FIG. 12. PDF of s for a Burgers run. Note the nearly flat slope at negative-s values. 
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